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\^ , Abstract 

O ■ 

(^ ■ Inspired by the multicanonical approach to simulations of first- 

C^^ , order phase transitions we propose for g-state Potts models a combi- 

nation of cluster updates with reweighting of the bond configurations 
Q^ ' in the Fortuin-Kastelein-Swendsen-Wang representation of this model. 

Numerical tests for the two-dimensional models with q = 7,10 and 
20 show that the autocorrelation times of this algorithm grow with 
p ' the system size V as t oc V°, where the exponent takes the optimal 

D . random walk value of q ~ 1. 
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1 Introduction 

Monte Carlo simulations of first-order phase transitions [jl[ in the canoni- 
cal ensemble are severely hampered by extremely large autocorrelation times 
r ex exp(2crL^~^) where a is the (reduced) interface tension between the co- 
existing phases and L^~^ is the cross-section of the system p|. To overcome 
this problem Berg and Neuhaus 0, |] have recently introduced multicanon- 
ical simulations which are based on reweighting ideas and can, in principle, 
be combined with any legitimate update algorithm. Using local update al- 
gorithms (Metropolis or heat-bath) it has been demonstrated in several ap- 
plications 1^ that the growth of autocorrelation times with system size is 
reduced to a power-law, r oc V^ with a > 1. For the two-dimensional g-state 
Potts model values of a ~ 1.3 have been reported for q* = 7 |P and g = 10 

i- 

Since by construction the multicanonical energy distribution is constant 

over the interesting energy range, invoking a random walk argument, one 
would expect an exponent a = 1 for an optimally designed update algo- 
rithm. The purpose of this note is to present for Potts models a cluster 
update variant of the multicanonical approach which is optimal in this sense. 
Basically the idea is to treat the cluster fiips in the first place and to reweight 
the bond degrees of freedom instead of the energy. 

2 The algorithm 

The basis of cluster update algorithms [0, H is the equivalence of the Potts 
model 

^Potts = E e"^""; ^ = - E ^'^.-^ ; cr, = 1, . . . , g, (1) 

with the Fortuin-Kastelein (FK) and random cluster (RC) representations 

-^Potts = ZpK = Z^c, (2) 

where 

^FK = E E n [p^a,<r,K„l + Kjfi] , (3) 

{o-i}{b,j} (ij) 



and 

with 

p = exp(/3) - 1. (5) 

Here hij = or 1 are bond occupation numbers and Nc{{hij}) denotes the 
number of connected clusters (including one-site clusters). According to (|^) 
a Swendsen-Wang cluster update sweep then consists in 1) setting hij = if 
(jj 7^ (Tj , or assigning values and 1 with relative probability 1 : p if a, = aj , 
2) identifying clusters of spins that are connected by "active" bonds {pij = 1), 
and 3) choosing a new random value 1 . . . g independently for each cluster. 

By differentiating In Z with respect to (3 it is easy to see that the average 
of the energy E = — J2{ij) ^aio-j can be expressed in terms of the average of 
the number of active bonds B = I](ij) %, 

^ = -(E) = i±l{B), (6) 

op p 

and for the specific heat per site C one finds 

,2 d{E) P+l,^, , fP+A\,r,2, /ON2 



CV/P^ = -^ = -^(5) + I^^J [{B^) - {By) . (7) 

Eq. (H) suggests that the bond histogram P^^^^{B) should develop for [3 = 
/3i ±5/5 a pronounced peak around Bo^d = — ^Eo^d, and for [3 ^ [3t a. double- 
peak structure similar to P^^l{E). In fact, as is illustrated in Fig. 1 for 
g = 7 and L = 60, a plot of P^^^'^ versus ^^5 is hardly distinguishable from 
P^^^{E). For other values of q and L the comparison looks very similar. 

In terms of P^-^n'^ the slowing down of canonical simulations is thus caused 
by the strongly suppressed configurations near the minimum between the two 
peaks, analogous to the well-known argument for P™^. To enhance these 
probabilities we therefore introduce in analogy to multicanonical simulations 
a "multibondic" partition function 

^mubo = H I] n [p^^^<yJK,^ + K,,q\ exp(-/bond(5)), (8) 

{o-»} {bij} (ij) 



where /bond(-B) = \nP^°j^'^{B) between the two peaks and /bond(-B) = oth- 
erwise. Of course, as in muhicanonical simulations, any reasonable approx- 
imation of P^°^'^{B) can be used in practice. Canonical expectation values 
can always be recovered exactly by applying the inverse of the reweighting 
factor exp(-/bond(5))- 

Obviously, once the bij are given, we can update the ai exactly as in the 
original Swendsen-Wang cluster algorithm. To update the bij we proceed 
as follows. If ai ^ aj then the bond bij is never active and we always set 
bij = 0. If (Jj = aj then we define B' = B — bij and choose new values 
bij = or 1 with probabilities Vibij = 0) = A/'exp(— /bond(-B')) and Vipij = 
1) = Afp exp(— /bond(-B' + 1)), where A/" is a trivial normalization factor 
(A/" = 1/ [exp(-/bond(5')) +p exp(-/bond(-B' + 1))] )• Since this is nothing 
but a local heat-bath algorithm for the bij the whole procedure is obviously 
a valid update algorithm. 

3 Results 

To evaluate the performance of the multibondic (mubo) cluster algorithm 
we performed simulations of the Potts model (|lD in two dimensions with 
g = 7, 10 and 20. The investigated lattice sizes and simulation temperatures 
are compiled in Table 1. For comparison we run for the same parameters 
also standard multicanonical (muca) simulations using the heat-bath update 
algorithm. In each run we recorded A^ = 100 000 measurements of E, B and 
two definitions of the magnetization in a time-series file. (The only exception 
is the multicanonical simulation for q = 7, L = 100 with A^ = 30 000.) 
Between the measurements we performed M lattice sweeps, with M adjusted 
in such a way that the autocorrelation times in units of measurements and 
thus the effective statistics of practically uncorrelated data was roughly the 
same in all simulations. 

To make sure that the new multibondic algorithm was implemented cor- 
rectly, we have analyzed some of the usually considered canonical quantities 
such as the specific heat C and the Binder parameter V = 1 — {E^)/3{E'^)'^. 
In Table 2 we compare results for the specific-heat maximum and Binder- 
parameter minimum for g = 10 obtained from our multicanonical and multi- 
bondic simulations. The error bars are estimated by the jackknife method 
[|T^. The data are in good agreement with each other and also with re- 



suits from independent canonical Metropolis |Tl| or single-cluster [0 high- 
statistics simulations. 

As discussed in Refs.|T^, Q in the context of a multicanonical multi- 
grid implementation it is not completely obvious which definition of the 
autocorrelation time should be used to characterize the dynamics of mul- 
ticanonical or multibondic simulations. One could, e.g., analyze the (mul- 
ticanonical/bondic) autocorrelation function of E, i?exp(/ene(-E)), B, or 
-Bexp(/bond(-B)), where fene{E) is the multicanonical analogue of /bond(-B)- 
More relevant from a practical point of view is the effective autocorrelation 
time [p!3| , |T^ for canonical observables which can be defined from the ratio 
of proper and naive error estimates (e/enaiv = v2t^) ■ The third possibility, 
which allows a direct comparison with previous work, is to define flipping 
(or more properly diffusion) times 4:T^^^ by counting the number of update 
sweeps that are needed to travel from E < -Emin to E > -Emax and back. Here 
-E'min.max (oi -Bmin.max for T^^^) are cuts which are usually chosen as the peak 
locations Eo^d{L) of the canonical probability distribution. Alternatively one 
could also use the infinite volume limits Eo^d of Eo^d{L) for all lattice sizes. 
For 2D Potts models this is straightforward since Eo^d are known exactly. In 
our simulations we have tested if E (or B) has passed the cuts after each 
sweep. We observed significantly larger t^'p when performing this test only 
every M'th sweeps, since then any cut-crossings during the M — 1 sweeps 
between the measurements cannot be detected. 

Our results for t^^ obtained in multicanonical and multibondic simula- 
tions for q = 7,10 and 20 are shown in the log-log plots of Figs. 2-4. Here 
we have used the canonical peak locations for the energy cuts. Let us first 
concentrate on the results for g = 7 in Fig. 2 where we have included for 
comparison the data from previous multicanonical simulations and also 
Rummukainen's results for his hybrid-like two-step algorithm which com- 
bines microcanonical cluster updates with a multicanonical demon refresh 
[ p!5[] . Both cluster update versions show qualitatively the same behavior and, 
for L > 20, perform much better than the standard multicanonical algorithm. 
From least-square fits to 

r^ip = aV (9) 

we estimate a ~ 1.3 for multicanonical heat-bath and a ~ 1 for multibondic 
cluster simulations; see Table 3 where we also give results for fixed energy 
cuts. Our results for t^^^ are almost indistinguishable from t^^^ which, recall- 



ing Fig. 1, is no surprise. Furthermore we have also measured the effective 
autocorrelation times and find that they are systematically smaller for both 
algorithms. 

Unfortunately, for g = 10 and 20 the situation is less favorable for the 
multibondic algorithm. While we still find an exponent of a ~ 1, the pref- 
actor in @ turns out to be so large that we can take advantage of this 
asymptotic improvement only for very large lattice sizes. As can be seen in 
Fig. 3, for g = 10 the cross-over happens around L = 50. Extrapolating to 
L = 100 we estimate that the multibondic algorithm would perform for this 
lattice size about 1.5 times faster than the standard multicanonical heat- 
bath. For q = 20 the same comparison clearly favors the standard algorithm 
for all reasonable lattice sizes - and we certainly cannot recommend the new 
algorithm for large q. 

4 Conclusions 

In summary, we have proposed for Potts models a combination of cluster 
update techniques with reweighting in the random bond representation and 
shown that this approach is feasible in practice. In fact, it is technically 
not more involved than the standard multicanonical approach and one lat- 
tice sweep takes about the same CPU time. Numerical tests for the two- 
dimensional g-state Potts model with g = 7, 10 and 20 show that the multi- 
bondic cluster algorithm is optimal in the sense that the exponent a in the 
power-law, t'^'p = aV^, is consistent with a = 1, the value one would expect 
in an idealized random walk picture. For g = 7 the multibondic algorithm 
clearly outperforms the standard multicanonical heat-bath algorithm. Com- 
pared with Rummukainen's hybrid-like two-step cluster variant the multi- 
bondic autocorrelation times are smaller for all lattice sizes by a roughly 
constant factor of 1.5. For larger values of g, however, the prefactor a turns 
out to be relatively large, rendering the new algorithm for reasonable lattice 
sizes more efficient than multicanonical simulations only for g < go with go 
somewhat above 10. 

The multibondic cluster algorithm may be of value for a wide range of 
investigations, since it can be applied to any systems where conventional 
cluster update techniques are applicable. 
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Tables 



Table 1: Simulation parameters: /Sq is the inverse simulation temperature, 



M„ 



and Mmubo denote the number of lattice sweeps between measure- 



ments, and |-E|min,max and B 
autocorrelation time t^^^. 



min,inax 



are the cuts used in the definition of the 



Q 


L 


Po 


-'"muca 


1 P 

l-C' min 


-C/ 1 max 


-'"mubo 


-Dmin 


-Dmax 




20 


1.284690 


10 


426 


644 


10 


310 


462 


7 


40 


1.291050 


50 


1801 


2542 


30 


1310 


1840 




60 


1.292283 


100 


4139 


5675 


70 


3003 


4118 




100 


1.293089 


100 


11692 


15672 


100 


8476 


11367 




12 


1.407380 


10 


120 


247 


20 


91 


186 




16 


1.415340 


20 


222 


439 


40 


166 


326 


10 


20 


1.418864 


30 


353 


676 


60 


270 


512 




26 


1.421642 


70 


614 


1137 


140 


467 


867 




34 


1.423380 


200 


1065 


1938 


400 


813 


1474 




50 


1.424752 


200 


2349 


4185 


400 


1797 


3182 




4 


1.577747 


5 


6 


32 


10 


5 


24 




6 


1.639809 


6 


17 


72 


12 


14 


56 




8 


1.665033 


12 


33 


124 


25 


28 


96 




10 


1.676647 


17 


56 


192 


35 


45 


151 


20 


12 


1.683517 


32 


79 


280 


65 


66 


215 




14 


1.688195 


67 


112 


357 


135 


92 


295 




16 


1.690278 


87 


154 


470 


175 


122 


384 




18 


1.692013 


125 


194 


593 


250 


155 


485 




20 


1.693698 


175 


231 


728 


350 


201 


595 



Table 2: Two-dimensional 10-state Potts model: Comparison of results 
for specific-heat maxima and Binder-parameter minima from multicanonical 
(muca) and multibondic (mubo) simulations. 



L 


alg. 


I3c 

1 ^max 




/^Vn,in 


T/ 
'^ mill 


12 
12 


muca 
mubo 


1.40621(28) 
1.40733(29) 


44.89(18) 
44.92(18) 


1.39256(29) 
1.39373(29) 


0.50379(75) 
0.50521(79) 


16 
16 


muca 
mubo 


1.41480(17) 
1.41451(15) 


73.68(25) 

73.87(27) 


1.40757(17) 
1.40729(15) 


0.52503(53) 
0.52447(61) 


20 
20 


muca 
mubo 


1.41837(11) 
1.418602(95) 


109.92(32) 
109.33(33) 


1.41392(11) 
1.414150(97) 


0.53541(44) 
0.53649(44) 


26 
26 


muca 
mubo 


1.421325(71) 
1.421506(54) 


177.61(44) 
177.91(43) 


1.418786(72) 
1.418969(54) 


0.54451(34) 
0.54453(34) 


34 
34 


muca 
mubo 


1.423313(36) 
1.423312(26) 


296.13(52) 
297.14(50) 


1.421873(36) 
1.421872(26) 


0.55001(24) 
0.54941(23) 


50 
50 


muca 
mubo 


1.424801(28) 
1.424834(21) 


627.79(97) 
627.6(1.1) 


1.424155(34) 
1.424188(21) 


0.55437(21) 
0.55439(22) 



Table 3: Results for the dynamical exponent a in multicanonical (muca) and 
multibondic (mubo) simulations from fits to t^^^ = aV", using L-dependent 
cuts defined by the canonical peak locations (ftmax) and fixed cuts at the 
infinite volume limits of Edo and Bdo (ctfix)- 





'^max 


ttfix 


q 


muca 


mubo 


muca 


mubo 


7 
10 
20 


1.27(2) 
1.32(2) 
1.26(1) 


0.92(2) 
1.05(1) 
1.09(1) 


1.53(2) 
1.43(1) 
1.46(1) 


1.02(2) 
1.12(1) 
1.18(1) 



Figure Headings 



Fig. 1: Canonical energy and bond histograms for q = 7, L = 60, and 

f3 = 1.292283. The bond histogram is plotted vs [{p + l)/p]B, where 
p = exp(/9) — 1. 

Fig. 2: Log- log plot of autocorrelation times r^''' of the energy vs lattice size 
for q = 7, using L-dependent energy cuts defined by the peak locations 
of the canonical energy distribution. 

Fig. 3: Same as Fig. 2 for g = 10. 

Fig. 4: Same as Fig. 2 for q = 20. 
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